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Abstract 

We introduce SDPB: an open-source, parallelized, arbitrary-precision semidefinite 
program solver, designed for the conformal bootstrap. SDPB significantly outperforms 
less specialized solvers and should enable many new computations. As an example 
application, we compute a new rigorous high-precision bound on operator dimensions 
in the 3d Ising CFT, Ao- = 0.518151(6), = 1.41264(6). 
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1 Introduction 


In [1], Rattazzi, Rychkov, Tonni, and Vichi showed how to bound Conformal Field Theory 
(CFT) observables using convex optimization. Shockingly, the resulting bounds are some¬ 
times saturated by actual CFTs, allowing high-precision computations of nonperturbative 
quantities [2-9]. This modern incarnation of the conformal bootstrap [10, 11] has been 
applied to numerous theories, and the list is growing [12-33]. 

With recent analytical work [34-40], we now understand in principle how to formulate 
many interesting bootstrap calculations, like studies of four-point functions of scalars, 
fermions, conserved currents, stress-tensors, and mixed correlators combining all of these 
ingredients. Such studies may shed light on the classihcation of critical points in condensed 
matter systems, the conformal windows of 4d gauge theories, the landscape of AdS string 
vacua, and more. Each study inevitably culminates in an optimization problem that only 
can be solved numerically at present. 

SDPB is a custom optimizer with three purposes: 

1. To enable the next generation of bootstrap studies involving conserved currents, stress 
tensors, and other complex ingredients. 

2. To help improve predictions for CFTs already isolated with the bootstrap, like the 3d 
Ising CFT, 3d 0{N) vector models, and others. 

3. To demonstrate the potential for dramatic improvement in current numerical boot¬ 
strap techniques and to encourage researchers in numerical optimization and computer 
science to contribute ideas and expertise. 

Custom optimizers have improved bootstrap calculations in the past. For example in [5], 
a custom-written linear program solver enabled a new high-precision calculation of critical 
exponents in the 3d Ising CFT, surpassing what was possible with out-of-the-box solvers 
like Mathematica [41] (used in the original work [1]), GLPK [42], and CPLEX [43].^ 

Unfortunately, linear programming is not applicable to systems of multiple correlators 
or operators with spin. These more complicated cases can be attacked with semidehnite 
programming [17, 9], for which previous studies have relied on the solvers SDPA [44, 45] 
and SDPA-GMP [46]. The study [9], in particular, pushed SDPA-GMP to its limits, with each 
optimization taking up to 2 weeks. By contrast, SDPB can perform the same optimization 
in 1-3 CPU-hours, or 4-12 minutes on a 16-core machine. 

A general bootstrap problem can be approximated as a particular type of semidehnite 
program called a “polynomial matrix program” (PMP). SDPB implements a variant of 
the well-known primal-dual interior point method for semidehnite programming [47-50], 
specialized for PMPs. Specialization and parallelization are SDPB’s advantages. We are 
optimistic that better designs and algorithms can be brought to bear in the future. 

^The solver in [5] is more accurately called a semi-infinite program solver. 
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In section 2, we describe PMPs and the design of SDPB. This discussion is mostly 
independent of the conformal bootstrap and should be comprehensible to readers without 
a physics background. 

In section 3, as an application of SDPB, we set a new world-record for precision of critical 
exponents in the 3d Ising CFT, using multiple correlators as in [9]. Readers interested solely 
in physics can skip to this section. We conclude in section 4. 

SDPB is open source and available online at https://github.coni/davidsd/sdpb. 


2 Design of SDPB 


2.1 Polynomial Matrix Programs 


SDPB solves the following type of problem, which we call a polynomial matrix program (PMP). 
Consider a collection of symmetric polynomial matrices 


m;(x) 


( P"n(a:) . 



■ ■ (a^) ) 


( 2 . 1 ) 


labeled by 0 < n < and 1 < j < J, where each element is a polynomial. Given 

h G M^, we would like to 


maximize b ■ y over y G 

such that Mj{x) + P 0 for all a: > 0 and I < j < J- 

The notation M P 0 means “M is positive semidehnite.” 

As we review in section 3, a wide class of optimization problems from the conformal 
bootstrap can be written in this form. Conveniently, PMPs can be translated into semideh¬ 
nite programs (SDPs) and solved using interior point methods. In the next few subsections, 
we perform this translation and describe an interior point algorithm for solving general 
SDPs. Subsequently, we make this algorithm more efficient by exploiting special structure 
in PMPs. 


2.2 Translating PMPs into SDPs 

Let us begin by translating the PMP (2.2) into a more standard semidehnite program of 
the following form: 

maximize Tr(Cy) + b ■ y over y G M^, Y G ^ 
such that Tr(A*F) -|- By = c, and (2.3) 

r P 0, 
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where 


c e 
B e 

Au...,Ap,C e S^. (2.4) 

Here, is the space of K x K symmetric real matrices, and Tr(H*H) denotes the vector 
(Tr(HiF),... ,Tr(HpF)) G The SDP (2.3) is similar to those treated by solvers like 
SDPA-GMP, except that it includes the variables y G called “free variables” because they 
are not constrained to be positive. The matrix C will eventually be set to zero, but we 
include it for generality. 

The hrst step is to relate positive semidehniteness of polynomial matrices to positive 
semidehniteness of a single matrix Y. Let qo{x),qi{x ),..., be a collection of polynomials 
with degrees 0,1,... (for example, monomials qm{x) = x^), and dehne the vector qsix) = 
{qo{x),... ,qs{x)). We call qm{x) a “bilinear basis” because products qm{,x)qn{x) span the 
space of polynomials. In particular, any polynomial P{x) of degree d can be written 

P{x) = Tt^si+i{YiQsi{x)) + xTt^s 2 +i{Y 2 Qs 2 {x)), (2.5) 


where 


Qs{x) = qs{x)qs{x)^, 

51 ^ L^/2J, 

52 = L(^-1)/2J, (2.6) 

and Yi,Y 2 are (underdetermined) symmetric matrices 

Yi G 5^1 

Y2 G 5^^+^ (2.7) 

For a symmetric mxm polynomial matrix M [x) of degree d, we can apply this construction 
to each element, 

M{x') Trjjij^+l ((Q^j^ (x) (S) Imxm)) T 3^Tr]]j52+l(b^(Q(52(^) ® Imxm)); 

Yi G 

e ^2.8) 

where now Yi acts on (g) and y 2 similarly acts on (g) ^^d each trace is 

over the hrst tensor factor. 

The translation of PMPs into SDPs relies on the following theorem: 

Theorem 2.1. M{x) is positive semidefinite for x G if and only if it can he written in 
the form (2.8) for some positive semidefinite Yi and Y 2 . 
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Proof. One direction is simple: choose a vector v G and consider the pairing v'^M{x)v = 
®vv'^))+ PT{Y 2 {xQi.^{x) ^vv"^)). Suppose Yi,h 2 are positive semidehnite and 
X > 0. Then since Qsi(x) <^vv'^ and xQs 2 (x) ®vv'^ are both positive semidehnite, it follows 
that v'^M(x)v > 0. The other direction is less trivial. It has been proven directly in [51] 
and also follows as a consequence of the “Biform Theorem” of [52], using the substitution 
X = and results of [53].^ □ 

Theorem 2.1 lets us rewrite our BMP constraints (2.2) in terms of a collection of positive 
semidehnite matrices Yi,, Y 2 J. We equate each polynomial matrix Mfi^x) to an expression 
of the form (2.8). Individual matrix elements can be isolated by taking the trace over 
M™'-’ with symmetrized unit matrices 






(2.9) 


This gives a set of polynomial equalities 


n 


Tr(r2,-i(Q5,,(a;) 0 E^^)) + Tt{Y2,{xQsJx) (g) P”*)), 


N 

dj = max deg(M"(a;)), 

n=0 

5ji = \_dj/2\, 

6,2 ^ L(^i-1)/2J. (2.10) 


Equality between polynomials of degree dj is equivalent to equality at dj + 1 points. Thus, 
evaluating (2.10) at points xq, ..., Xd^, we obtain a set of affine relations between the yn and 
Y ^ 

Pj,rs{Xk) + ynP;',rs{Xk) = Tr (^2,-1 (Q^.i (Xfc) 0 P”^)) + Tr(F2, (Xfc) ® P”*)) 

n 

Q < r < s < rrij, /c = 0,..., dj. (2-11) 


Let us group the 1^’s into a single block-diagonal positive semidehnite matrix 


Y = 


(Yi 0 
0 ^2 


\0 0 


The constraints (2.11) now take the form 


0 \ 
0 

Y2JJ 


( 2 . 12 ) 


TT:{ApY) + {By)p = Cp. ( 2 . 13 ) 

^We thank Pablo Parrilo for pointing this out. 

^We could alternatively match coefficients on both sides of (2.10) as in [17, 7, 9]. We will see in 
subsection 2.5.2 why pointwise evaluation is preferable. 
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where p runs over all tuples {j,r,s,k) satisfying 0 < r < s < rrij, and 0 < k < dj. The 
matrices Ap,B,C and vector c are given by 

/O ■■■ 0 0 ■■■ 0\ 


A 


(j,r,s,k) 


0 0 ■■■0 

0 ■■■ 0 XkQsj2i.^k) ® ■■■ 0 ’ 


(2.14) 


\0 ■■■ 0 

B(j^r,s,k),n = ~Pj,rs{.^k)i 

C(j,r,s,k) — Ej,rs{^k)} 

C = 0. 


0 ■■■ 0 / 

(2.15) 

(2.16) 

(2.17) 


This completes the translation of our PMP into an SDP (2.3). Note that the matrices 
Ap are far from generic. Exploiting this fact will help us solve PMPs much more efficiently 
than a generic SDP. SDPB is specihcally designed to solve SDPs with constraint matrices of 
the form 


Aj,r,s,k) = ^b{Vi^kvlk®E''^), (2.18) 

b € blocks^ 

where Bb{M) denotes the block-diagonal matrix with M in the &-th block and zeros every¬ 
where else, 


/o 


Bb{M) = 


0 ■■■ 0 \ 


0 M 0 


Vo 


0 


(2.19) 


0 / 


and the sets blocks^- are disjoint for different j. In the example above, we have 

blocks^- = {2j-l,2j}, 

—1,/c 

V2j,k = xl^^qs^^ixk). 


( 2 . 20 ) 


2.3 Semidefinite Program Duality 

Duality plays an important role in our interior point algorithm, so let us briefly review it. 
The problem (2.3) is related by duality to the following “primal” optimization problem: 

V : minimize c • x over x G X G , 

such that X = X]p=i ApXp - C, , . 

B^x = b, ^ ’ 

X yo. 
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We refer to the problem (2.21) as V for “primal” and (2.3) as V for “dual.” We say that V is 
“feasible” if there exist x, X satisfying the constraints (2.21). Similarly, V is feasible if there 
exist y,Y satisfying the constraints (2.3). The “duality gap” is dehned as the difference in 
primal and dual objective functions, c ■ x — Tr(CF) — b ■ y. 

For our purposes, the statement of duality is as follows: 

Theorem 2.2 (Semidehnite Program Duality). Given a feasible point {x,X) of V and a 
feasible point {y,Y) ofV, the duality gap is nonnegative. If the duality gap vanishes, then 
{x,X) and {y,Y) are each optimal solutions ofV and V, and furthermore XY = 0. 


Proof. Suppose we have feasible solutions {x,X) and {y,Y). The duality gap is given by 


c - X — Tr(Cy) -b-y = c - x — Tt 


^ApXp-X^ yj -b-y 


= c - X + Tr(Xy) — X ■ Tr(74*y) —b-y 
= C-X + Tr(Xy) - X • (Tr(A,y) + By) 
= Tt{XY) > 0, 


( 2 . 22 ) 


where nonnegativity follows because X and Y are positive semidehnite. Now suppose 
Tr(Wy) vanishes. Clearly this implies XY = 0 identically (this condition is called “com¬ 
plementarity”). Because c • x is bounded from below by the dual objective Tr(Cy) + b - y 
and also equal to the dual objective, the point (x, X) must be optimal. Similarly, since 
Tr(Cy) + b - y is bounded from above by the primal objective c • x and also equal to the 
primal objective, the point {y, Y) must be optimal as well. □ 


Unlike in linear programming, there is no guarantee that either V or T> will attain 
their respective optima, or that the duality gap will vanish. For this, we need additional 
regularity assumptions. One of them is Slater’s condition, which says that the duality gap 
vanishes if there exist strictly feasible solutions to the primal and dual constraints — i.e. 
solutions where X, U >- 0 are positive-dehnite. Slater’s condition is generic in the sense that 
a small perturbation of a feasible but not strictly-feasible problem will typically produce a 
strictly-feasible problem. 


2.4 An Interior Point Method 

The idea behind primal-dual interior point methods is to solve both the primal and dual 
equations simultaneously to hnd an optimal point q = {x, X,y,Y). As we saw in the 
previous subsection, the optimum is (generically) achieved by a pair of positive semidehnite 
matrices X, Y satisfying the “complementarity” condition XY = 0. Most algorithms work 
by deforming this condition to 


XY = yl 


(2.23) 


for some nonzero /r, where I is the identity matrix. The constraints together with (2.23) then 
have a nniqne family of solntions called the “central path:” g(/i) = (x(/i), X(/i), |/(p), y(/i)) 
indexed by p G By following the central path from positive p towards /i = 0, we can 
hnd the optimnm of the original problem. 

In practice, instead of moving precisely along the central path, we nse the following 
strategy. Consider an initial point q = {x, X, y, Y) snch that X, Y are positive semidehnite. 
Our goal is to decrease Ti{XY) and move towards the constraint locus while maintaining 
positive semidehniteness. 

• Set p = l3Ti{XY)/K for some /9 < 1, where K is the number of rows of X. 

• Use Newton’s method to compute a direction dq = {dx, dX, dy, dY) towards the 
central path with the given y. 

• Take a step along dq, taking care not to violate the positive semidehniteness of X, Y. 
This should result in a reduction of Tr(XU) by roughly a factor of 13. 

• Repeat. 

This is essentially Newton’s method with a moving target. 

An important advantage of this method is that the initial starting point [x, X, y, Y) 
need not satisfy any of the equality constraints in (2.3) and (2.21). As long as we start with 
positive semidehnite X, Y, and the problem is feasible, the above method will converge to 
a point where the equality constraints are satished. 


2.4.1 Newton Search Direction 

Let us describe a single Newton step in more detail. The direction dq is dehned by replacing 
q ^ q + dq and solving the constraint equations (2.3, 2.21) and complementarity equation 
(2.23) at linear order in dq, 


X + dX 

B^{x + dx) 
Tt{A,(Y + dY)) + B{y + dy) 
XY + XdY + dXY 


'^Ai^Xi + dxi) - C, 

i 

b, 

c, 

fil. 


The residues 


P = Y,AiXi-X-C, 

i 

= b — B^x, 

= c-AT{A,Y)-By, 
= yl-XY, 


(2.24) 


P 

d 

R 
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(2.25) 


measure the failure of the current point to satisfy the constraints. These residues will 
decrease with each Newton step. The linearized equations (2.24) can then be written 



^-d - Tt{A,Z)^ 
P + y^Aidxi, 

i 

X-\R-dXY), 


where Z = X ^{PY — R) and the “Schur complement” matrix S is dehned by 

Sij = Ai{AX-^AjY). 


(2.26) 

(2.27) 

(2.28) 


(2.29) 


To hnd the search direction, we hrst solve (2.26) for dx, dy, and then plug into (2.27) and 
(2.28) to determine dX, dY. Naively applying (2.28) leads to a dY that is not necessarily 
symmetric, which would take us outside of the domain of dehnition of Y. Several solutions 
to this problem have been proposed [54]. Our approach, following [49, 55, 56], will be to 
symmetrize dY by hand, replacing (2.28) with 

= X-\R-dXY), 

dY = + dY^) ■ (2.30) 


2.4.2 Mehrotra Predictor-Corrector Trick 

The most expensive operations in the search direction calculation are forming the Schur 
complement matrix S and solving the Schur complement equation (2.26). We’d like to 
perform them as rarely as possible. A simple modihcation to the naive Newton’s method, 
due to Mehrotra [57], allows us to get closer to the central path while reusing S and most 
of the work done in solving (2.26). 

The rough idea is to solve the constraint and complementarity equations at higher 
order. This proceeds in two steps, called the “predictor” and “corrector” steps, respectively. 
For the predictor step, we compute a direction as described above, which we call dqp = 
{dXp,dXp,dyp,dYp). We then replace the linearized complementarity equation (2.24) with 

XY + XdY + dXY + dXpdYp = yl (2.31) 

and re-solve to obtain a corrector direction dqc- In the corrector step, we may (optionally) 
use a smaller deformation parameter /r —)■ /Xc to get closer to Tt{XY) = 0. Note that 
the replacement of (2.24) with (2.31) does not change the Schur complement matrix Sij, 
so it can be reused, together with any matrix decompositions performed in solving (2.26). 
Altogether, the corrector step amonnts to simply replacing 


R^ fij -XY - dXpdYp 


(2.32) 


before solving (2.26), (2.27), and (2.30). 


10 


2.4.3 Termination Conditions 


We say a point q is “primal feasible” if the residues p, P are sufficiently small. Similarly, the 
solution is “dual feasible” if the residue d is sufficiently small. The precise conditions are 

primal feasible: primalError = maxjj{|pj|, \Pij\} < primalErrorThreshold; 
dual feasible: dualError = maxi{|(ij|} < dualErrorThreshold, 

(2.33) 

where primalErrorThreshold <C 1 and dualErrorThreshold 1 are parameters chosen 
by the user. 

An optimal point should be both primal and dual feasible, and have (nearly) equal 
primal and dual objective values. Specifically, let us define dualityGap as the normalized 
difference between the primal and dual objective functions 

IprimalObjective — dualObjective| 
max{l, IprimalObjective + dualObjective|}’ 

c • X, 

Tr(Cy) + b-y. (2.34) 

if 

dualityGap < dualityGapThreshold, (2.35) 

where dualityGapThreshold 1 is chosen by the user. 


dualityGap = 

primalObjective = 
dualObjective = 

A point is considered “optimal” 


2.4.4 Complete Algorithm 

Our complete interior point algorithm is as follows. 


1. Choose an initial point q = {x, X,y,Y) = (0, OpJ, 0, Ox^-f) where flp and flx> are real 
and positive. This point probably does not satisfy the constraints. 

2. Compute the residues (2.25) and terminate if q is simultaneously primal feasible, dual 
feasible, and optimal. (Sometimes we may wish to use a different termination criterion, 
see below.) 

3. Let fi = Tt{XY)/K and compute the predictor deformation parameter pp = (5py 
where 





Here, Anfeasibie € (0,1) is chosen 


if q is primal and dual feasible; 
otherwise. 

by the user. 


(2.36) 
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4. Compute the predictor search direction dqp by solving eqs. (2.26, 2.27, 2.30) with 
R = IJpI - XY. 


5. Compnte the corrector deformation parameter jUc = as follows. Let r = Tr((X + 
dXp){Y + dYp))/{iiK) and set 


/S 

Pc 


if r < 1; 
r otherwise. 


min(max(/3feasibie,/3), 1) 
max(/3infeasible) /3) 


if q is primal and dnal feasible; 
otherwise. 


(2.37) 

(2.38) 


where /^feasible £ (0,1) is a parameter chosen by the nser. This choice of Pc is modeled 
after the one in SDPA. 

6. Compnte the corrector search direction dqc by solving eqs. (2.26, 2.27, 2.30) with 

R = lj-XY-dX,dY,. 

7. Determine the primal and dnal step lengths 

a-p = min(7a(X, dXc), 1), 

av = min(7Q!(F, (717), 1), (2.39) 

where a{M,dM) is the largest positive real nnmber snch that M + a{M, dM)dM is 
positive semidehnite, and 7 G (0,1) is a parameter chosen by the user. 


8. Replace 

X ^ X + a-pdxc, 

X —^ X T Oi'pdXc, 
y ^ y + apdyc, 

Y Y + avdYc, (2.40) 

and go to step 2. Note that the replacement (2.40) is guaranteed to preserve positive 
semidehniteness of X and Y. 


If the current point is close enough to a primal (or dual) feasible region, the step-length 
ap {ap) in (2.39) can be exactly 1. When this occurs, the replacement (2.40) will exactly 
solve the primal (dual) equality constraints, up to numerical errors. This follows from 
linearity of the equality constraints, together with the fact that symmetrizing Y in (2.30) 
has no effect on the constraints. In cases where we only care about primal or dual feasibility, 
the iteration can be stopped here, see section 3.4. 


2.5 Specialization to Polynomial Matrix Programs 

As mentioned in section 2.4.2, the most expensive part of the search direction calculation 
is computing the Schur complement matrix Spg = Ti{ApX~^AqY) and inverting (2.26) to 
obtain dx, dy. In this section, we will specialize to PMPs, and study how these calculations 
can be made more efficient. For similar optimizations, see [58]. 
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2.5.1 Block Structure of the Schur Complement Matrix 


Recall that for PMPs, the matrices Ap are given by (2.18) with the index p running over 
tuples {j,r,s,k) satisfying 0 < r < s < and 0 < k < dj. Since X and Y have the 
block structure (2.12), the Schur complement matrix Sp^p^ is block-diagonal: it has nonzero 
entries only if ji = j 2 , 


S = 


0 

0 5 ( 2 ) 


0 \ 
0 


(2.41) 


\ 0 0 5(‘^)/ 

The size of the j-th block dim5(-^) is equal to the number of choices for (r, s, fc), 

dims® = + 1 ). 


(2.42) 


Now consider equation (2.26), 

T ~ Tr(A*Z)^ 


where 




(2.43) 


We could solve it using an LU (lower triangular x upper triangular) decomposition of T, 
but this is extremely expensive, taking cubic time in diniT = dim5(-^) -|- N. 

We should use the block structure of T to our advantage. Let S = LL^ be a Cholesky 
decomposition of S (which can be computed efficiently for each block 5(-^) = L('^)L('^)^), and 
consider the decomposition 



f L 0\ fl 0 \ 

1J 1^0 B^L-^L-^B y) 0 1 )■ 


(2.44) 


The outer matrices on the right hand side are triangular, and can be solved efficiently by 
forward/backward-substitution. Meanwhile, the matrix Q = B^L~^L~^B typically has a 
much smaller dimension than S, dimQ = N dimT, so the middle block-matrix can be 
easily solved using a Cholesky decomposition.^ 

Unfortunately, the decomposition (2.44) is numerically unstable when S is ill-conditioned. 
Indeed, suppose S has a very small eigenvalue (so L does too), while the full matrix T does 
not. Then quantities like L~^B that appear in (2.44) will have large entries which must 
nearly cancel when recombined into a solution of (2.43). Near-cancellation reduces numerical 
precision. 

These problems stem from the off-diagonal blocks of T, which ultimately come from the 
free variables y in our semidehnite program. Several authors have considered the problem 

"‘The matrix Q, is often called the “Schur complement” in the block matrix decomposition (2.44). We 
will continue to use the words “Schur complement” to refer to S (which is the Schur complement of a 
different block matrix system). Hopefully this will not cause confusion. 
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of efficiently and stably solving semidefinite programs with free variables, with no obvious 
consensus [59]. For example, [60] suggests eliminating free variables by explicitly solving 
the primal constraint B^x = b and taking appropriate linear combinations of the matrices 
Ap. However this procedure destroys the sparsity structure of S, making it no longer block 
diagonal and forcing us to use an expensive full Cholesky decomposition. 

Preserving the structure of S and T is paramount. The simplest way to stabilize (2.44) 
is to increase the precision of the underlying arithmetic. In practice, this appears to 
be necessary anyway for larger-scale bootstrap problems, see appendix A. For additional 
stabilization, we use an old trick of adding low-rank pieces to S to make it better-conditioned. 
Suppose Ui are vectors, each with a single nonzero entry, such that S' = S + ^ • Uiuf = 
S -|- UU'^ has no small eigenvalues.^ Note that S' differs from S in only a few diagonal 
entries — in particular it has the same block structure. Now let us replace (2.43) with the 
system 


/dx\ /S' -B -U\ /dx\ /-d-TT{A,Z)\ 
r I d|/ j = ( 0 0 j ( di/ j = I p j . (2.45) 


By solving for 2 ; and substituting back in, it is easy to see that (2.45) is precisely equivalent 
to (2.43). However, the advantage is that because S' is well-conditioned, a decomposition 
like (2.44) is numerically stable. Indeed, dehning B' = {B U), we have 


T' 


Q' 


f L' 0\( 
B'^L'-^L'-^B' - 


1 0 
0 Q' 

'0 0 
.0 1 


0 





(2.46) 

(2.47) 


where S' = L'L'^. Because Q' is no longer necessarily positive semidehnite, we are forced 
to use an LU decomposition to invert the middle matrix in (2.46), which is slightly more 
expensive than Cholesky decomposition. Fortunately, Q' is usually small enough that this 
cost is inconsequential. If T itself is ill-conditioned, then this will manifest as instabilities 
when we try to LU decompose Q'. In this situation, there is little we can do to avoid 
imprecision. 


2.5.2 Computing the Schur Complement Matrix 

Now that we know what to do with the Schur complement matrix S', let us compute it 
efficiently. The fact that Ap has low rank is helpful. This explains why we chose to evaluate 
the polynomial equalities (2.10) at discrete points Xk in (2.11). Matching polynomial 
coefficients on each side, as done in [17, 7, 9], leads to higher-rank matrices Ap. The 

^The Ui can be found as follows. During Cholesky decomposition S = , we keep track of the 

geometric mean A of the diagonal entries. Whenever we encounter a diagonal entry La < 6A, where 9 1 

is a parameter, we replace La —>■ La + A, which amounts to choosing Ui = Aci where is a unit vector in 
the i-th direction. 
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helpfulness of low-rank constraints in solving SDPs, and their appearance in polynomial 
optimization, is well known [58]. 

Recall that X and Y are block diagonal (2.12), with each block Xf, acting on a tensor 
product of the form (gjR™. Let e ]ij('5+i)x(<5+i) (denote the (r, s)-th block of X^ in 

the second tensor factor, which acts on 


Since S is block diagonal, we need only compute elements with ji = j 2 = j- We have 


% 


,ri,si,ki),{j,r2,S2,k2) 


b € blocksj 


7 . 

b € blocksj 




+ (ri -H- Si) + (r2 -H- S 2 ) + (ri -H- Si, rs -H- S 2 )) -(2.48) 

Thus, instead of performing repeated matrix multiplications to calculate Tt{Ap^X~^A p^Y), 
we can precompute the bilinear pairings 


U, 


(b) 

{dj +l)s+/i:i ,(dj+l)r -\-k2 




ib) 


{dj+l)s+k2,{dj+l)r+ki 




(2.49) 

(2.50) 


and plug them into (2.48). Whereas forming S is often the most expensive operation in less- 
specialized solvers, the method outlined here makes it subdominant to other computations, 
see appendix B.l. 


2.5.3 Computing Step Lengths 

To hnd step lengths a-p, ap,, we must be able to compute a{M, dM) where M is a positive 
semidehnite matrix and a(M, dM) is the largest positive real number such that M -|- 
a{M, dM)dM is positive semidehnite. Let M = LL^ be a Cholesky decomposition of 
M. Since M + adM = L{1 + aL~^dML~'^)L'^, we have 

a{M,dM) = —l/min-eigenvalue(L“^(iML“^). (2.51) 

In SDPB, we compute all the eigenvalues of L~^dML~"’" using a QR decomposition and then 
simply take the minimum. Some solvers, like SDPA, implement the more efficient Lanczos 
method [61] for computing the minimum eigenvalue. In practice, the step-length calculation 
is a small part of the total running time, so we have neglected this optimization. 


2.6 Implementation 

SDPB is approximately 3500 lines of C-I--I-. It uses the GNU Multiprecision Library (GMP) 
[62] for arbitrary precision arithmetic, and MPAGK [63] for arbitrary precision linear alge¬ 
bra. The relevant MPAGK source hies are included with SDPB, with some slight modihca- 
tions: 
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• The Cholesky decomposition routine Rpotrf has been modihed to implement the 
stabilization procedure described in footnote 5. 

• Some loops in the LU decomposition routine Rgetrf have been parallelized. 

SDPB also depends on the Boost C++ libraries [64] and the parsing library tinyxml2 [65]. 

Previous experience shows that high-precision arithmetic is important for accurately 
solving bootstrap optimization problems. It is not fully understood why. The naive reason 
is that derivatives d^d^gi\,e{z,z) of conformal blocks vary by many orders of magnitude 
relative to each other as A varies. It is not possible to scale away this large variation, 
and answers may depend on near cancellation of large numbers. In practice, the matrix 
manipulations in our interior point algorithm “leak” precision, so that the search direction 
{dx, dX, dy, dY) is less precise than the initial point (x. A, ?/, Y). By increasing the precision 
of the underlying arithmetic, the search direction can be made reliable again. This strategy 
(which we adopt) comes at a cost of increased runtime and memory usage. Better strategies 
for dealing with numerical instabilities in bootstrap problems could bring enormous gains. 

SDPB is parallelized with OpenMP [66]. Because most matrices appearing in the interior 
point algorithm are block-diagonal, most computations are “embarrassingly parallel:” dif¬ 
ferent blocks can be distributed to different threads. (The most prominent exception is the 
LU decomposition of Q', which is why Rgetrf was modihed.) Consequently, performance 
scales nearly linearly with the number of cores, as long as the number of matrix blocks 
is sufficiently large. This is usually the case for interesting bootstrap problems, where J 
(which sets the number of blocks) is typically much larger than the number of available 
cores. It should be possible to achieve favorable scaling up to dozens or even hundreds of 
cores using MPI and more careful memory management. Further scaling should be possible 
with more hne-grained parallelization. 

SDPB is available online at https://github.com/davidsd/sdpb under the MIT license. 
The source code is carefully commented and written for readability (to the extent that C++ 
code is ever readable). We hope this will encourage customization and improvement. 

3 Example Application: 3d Ising Critical Exponents 

3.1 A 3d Ising Optimization Problem 

Bootstrap optimization problems can be naturally approximated as PMPs [17, 9]. In this 
section, we review the PMP for the system of correlators {(acrcra), (acree), (eeee)} in the 3d 
Ising CFT. We will be brief. Much more detail is given in [9]. 

Associativity of the Operator Product Expansion (OPE) for {(crcraa), (acree), (eeee)} 
implies the consistency condition 


(1 



(3.1) 
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Here, 0~^ runs over Z 2 -even operators of even spin and 0~ runs over Z 2 -odd operators of 
any spin. We have separated out the unit operator. A and £ are the dimension and spin of 
O, respectively. The object V-^a,£ is a 5-vector and is a 5-vector of 2 x 2 matrices 




0 

0 




V-i 


-i)~4iXAn,v) 


u,v) 

Un. 

-lYFZZMj 


y+A/ = 


0 0 
0 0 
0 fyxMv) 
0 0 
0 0 


\ 


( ° 

kFlXei^,v) 

\\FlTMv) 

0 

( ° 

lFZZ{u,v) 


0 


(3.2) 


with entries that are functions of conformal cross-ratios u and n, 


'M 






= V 


^ijAkl , 

9a,£ ( 


± u 


A ■ ■ A. 

9A,i 


‘^ij — 


A, - Aj. 


(3.3) 


The conformal blocks, which are known special functions. 

The OPE coefficients XaaO, Xaeo, X^eo and dimensions A are not known a priori. Nonethe¬ 
less, we can constrain them by understanding when it is possible for the terms in (3.1) to 
sum to zero. To do this, consider a 5-vector of functionals a = {a ^,..., a^), where each a* 
acts on the space of functions of u and v. Acting on (3.1) with a gives 


(1 1) 3. y,„,„(()+ 5; (A 

A / Q + 


Wo-O -^eeo) <3 • y+A/ 


X(j(tO 

XeeO 


• y-A,f- = 0- 


o- 


(3.4) 


The bootstrap logic, in the spirit of [1], is as follows. First we make an assumption about 
which A,£ appear in (3.4). We then search for a functional a such that 


(1 


uao (;) 

> 

0, 



•3 • y+A,(- 


0, 

for all Z 2 -even operators with even spin. 


<3 • V-A/ 

> 

0, 

for all Z 2 -odd operators of any spin. 

(3.6) 


The OPE coefficients Xijk are real in a unitary GET. Thus, if a exists, then it is impossible 
to satisfy the consistency condition (3.1), and our assumption is ruled out. By making 
different assumptions and searching for functionals a, we can map out the space of allowed 

A. 
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3.2 Approximation as a PMP 


The conditions (3.5) define a feasibility problem with an infinite number of semidefiniteness 
constraints (one for each A,i). To obtain a PMP, we choose a particular type of functional 

ai/]= E (3-6) 

m>n 

m-\-n<A 


where u = zz, v = (1 — z){l — z). Although the bootstrap logic does not depend on the 
types of functionals considered, only functionals of the form (3.6) lead to a PMP. Other 
types of functionals require different optimization methods. 

Derivatives of conformal blocks have a systematic approximation in terms of positive 
functions times polynomials. 


am an Ai2,A34 

%i/Ay 


{z,z) 



Xe{A)p 


Ai2,A34;mn 

I 


(A), 


(3.7) 


where are polynomials and ^^(A) are functions that are positive for all A in 

a unitary CFT. It follows that 




ij,kl 


z-^ ib,A,£V 




(3.8) 


where P^^^^’”^”(A) are linear combinations of p^‘-”^''*’”*”'(A). Using this approximation, and 
stripping off the positive factors ^^(A), (3.5) becomes a PMP: 


find such that: 


(1 1)Z„(0) 


> 0 , 


ZiiA) y 0, for all Z 2 -even operators with even spin, 
Yi{A) > 0, for all Z 2 -odd operators in the spectrum. 


Here ^^(A) are polynomials and Ze{A) are polynomial matrices in A defined as 


>^KA) 

Ze{A) 


mn 


P' 


cxe.aeimn / 


(A) + 


i Tjea,ae]mn 


(A)], 


El 


a 


. . .(j(j,crcr;mn / a \ 

(A) 

4 Tyaa,€e;mn/ \\ . 5 7-)crcr,ee;mn 


r/’ P 
2 — 


‘w+<nPi:r' 


(A)) 


b 


(3.9) 




ai„P-J;r’""(A} 


( 3 . 


Typically, we assume that A can vary arbitrarily above some minimum value Ajnin(^)- 
Writing A = Amin(-^) + x, we have positive semidefiniteness for all a: > 0. 

There are two important differences between (3.9) and our original PMP (2.2): 
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1. In (2.2) we have a finite number of positive semidefiniteness conditions j = 1,..., J, 
whereas here we have an infinite number since i can be any nonnegative integer. In 
practice, we include spins i up to some large but finite £max- As long as £max is large 
enough, a functional obtained by solving the problem with i < £max should also satisfy 
positive semidefiniteness for spins i > f'max- The proper choice of £ max depends on the 
problem at hand, see appendix A. See [27] for a more careful analysis. 

2. In (2.2) we are trying to optimize an objective function b ■ y, whereas here we are 
only interested in feasibility. To determine feasibility, we can pick the trivial objective 
function b = 0 and run our interior point algorithm until it becomes dual feasible. 


3.3 Setting Up SDPB 


The natural objects entering our calculation are (approximate) derivatives of conformal 
blocks Xf(A)pf^^’^®^’™'"(A), as opposed to just the polynomials pf^^’^®'‘’™'"'(A). Removing 
positive factors does not affect positive semidefiniteness, but it does affect the scaling of the 
resulting SDP. Restoring quantities to their “natural” size can improve numerical stability 
and performance.® SDPB provides a few different ways to implement this rescaling. 

As we saw in section 2.2, translating a PMP into an SDP requires a bilinear basis qm{x). 
SDPB allows a choice of bilinear basis qm\x) for each j = 1,..., J. For bootstrap problems, 
we take qm\x) to be orthogonal polynomials with respect to the norm 

poo 

= / dxxei^mmii)+x)p{x)q{x), (3.11) 

Jo 

where i is the spin corresponding to j. The change of basis between orthogonal polynomials 
with respect to (•, ami monomials x™ (used in previous bootstrap applications of SDP) 
is extremely ill-conditioned at high degree. So although the choice of qm is unimportant in 
principle, it can have a dramatic effect on numerical stability.^ 

SDPB also requires a set of sample points at which to evaluate polynomials, as well 
as scaling factors that modify the constraints (2.15, 2.16) as follows: 


B 


{j,r,s^k),n 


C(j,r,s,k) 



(3.12) 


Additionally, the A(^j^r,s,k) are given by (2.18) with 


V2,.t = (3.13) 


similar observation was made for the algorithm in [5]. 

^We thank Pablo Parrilo for pointing out the usefulness of orthogonal polynomials in improving the 
numerical stability of polynomial optimization. 
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This Sfc-dependent rescaling gives an isomorphic, but potentially more numerically stable 
SDP. For bootstrap problems, it is natural to pick where i corresponds to j. 

The can be any sequence of distinct points. A natural choice are zeros of one of the qm 
of sufficiently high degree. 

To summarize, SDPB depends on the following input: 

• for each j = 1,..., J: 

— polynomial matrices M°(x),..., M^{x) of maximum degree dj, 

— bilinear bases qm\x) (m = 0,..., [(ij/2j), 

— sample points x^^'^ (/c = 0,..., dj), 

— sample scalings (A: = 0,..., dj), 

• an objective function b G 

SDPB reads this data in an XML format described in the manual. A Mathematica package 
that translates the above data from Mathematica expressions into XML is included with the 
source distribution. An example 2d bootstrap computation is also included. More details 
about SDPB’s input and output formats and its various settings can be found in the manual. 

3.4 Results 

As an application of SDPB, let us improve upon the determinations of 3d Ising critical 
exponents in [5, 9]. We make an exclusion plot for the operator dimensions (Ao-,A£) as 
follows. Fix (Ao-, Ag) and use SDPB to determine if the PMP (3.9) is dual-feasible, i.e. 
whether there exist {y, Y) satisfying their associated constraints. If the PMP is dual-feasible, 
then the given (Aq-, A^) are excluded. If the PMP is not dual-feasible, then we cannot 
conclude anything about (Ao-,Ae). By scanning over different points, we map out the 
excluded region in (Ag-, AJ space. 

To determine dual feasibility, we use a vanishing objective function 6 = 0 and run SDPB 
with the option —f indDualFeasible. This terminates the solver if {y,Y) are found satis¬ 
fying their constraints to sufficient precision (i.e. if dualError < dualErrorThreshold).® 
In practice, if SDPB finds a primal feasible solution {x,X) after some number of iterations, 
then it will never eventually hnd a dual feasible one. Thus, we additionally include the 
option —findPrimalFeasible, which terminates the solver whenever primalError < 
primalErrorXhreshold. The termination status of SDPB then determines whether a point 
is allowed or not: 

found dual feasible solution (A^., A^) disallowed, 

found primal feasible solution (Ao-,Ag) allowed. ^ ' 

The precise SDPB options used for the computations in this work are described in appendix A. 

®If we kept running the solver, primalObjective would converge towards dualObjective = 0 and an 
optimum would eventually be reached. 
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A, 


allowed region for A = 19, 27, 35,43 



Figure 1: Allowed region for a Z 2 -symmetric 3d CFT with two relevant scalars, computed 
using SDPB with the system of correlators (aaaa), (aaee), and (eeee). The blue regions 
correspond to A = 19,27,35,43, in decreasing order of size. The larger black rectangle shows 
the current most precise Monte Carlo determinations of critical exponents in the 3d Ising 
CFT [67]. The smaller black rectangle shows the estimate for (Ag-, Ag) using c-minimization 
at A = 41 for the single correlator {aaaa) [5] . 


In figure 1, we plot the allowed regions for different numbers of derivatives labeled by 
A = 19,27,35,43,® corresponding to functionals a of dimension 275, 525, 855, and 1265, 
respectively.^® We focus on (Ao-, A^) near the 3d Ising CFT, leaving wider exploration to 
the future. The allowed region is an island that shrinks rapidly with increasing A.^^ The 
largest island, corresponding to A = 19 is the same as the allowed region in figure 5 of [9]. 
We can estimate the point towards which the islands shrink as follows. Let (oa, ^a) be the 
bottom-left point of the A-allowed island, and similarly let (ca, d^) be the top-right point. 
Dehne 


E^{r) = stddevAg{i9,27,35,43}(’^aA + (1 - r)cA), 

Ey{r) = stddevAe{i9,27,35,43}(?^^A + (1 - ?^)dA), (3.15) 

®Umax = 10,14,18, 22 in the notation of [9]. 

performance analysis for different values of A is given in appendix B.2. 

^^Each allowed region plotted in this work was computed by testing a grid of points and fitting curves to 
the boundary between allowed and disallowed gridpoints. The raw gridpoint data is available on request. 
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and let r^, Vy be the minima of E^, Ey respectively. Our estimate is^^ 

(A^, Ae) (r^a43 + (1 - r^)c43, ryb^^ + (1 - ry)d43) 

^ (0.5181478(5), 1-412617(4)), (3.16) 

where the errors are given by Ex{rx), Ey{ry). 

The dimensions (Ao-,Ag) were estimated in [5] using the conjecture that the 3d Ising 
CFT minimizes c = A^/Ao-o-t^^ subject to the constraints of unitarity and crossing symmetry 
of (craacr). This conjecture, called “c-minimization,” is expected to be equivalent to the 
assumption that the 3d Ising CFT lives precisely at the kink in the dimension bound 
coming from the single correlator {aaaa). Although unproven, c-minimization’s advantage 
is that it allows one to estimate (A^., AJ using a single scalar correlator {aaaa). Bootstrap 
computations for a single scalar correlator can be made relatively efficient with a modihed 
primal simplex algorithm [5, 70]. 

An advantage of multiple correlators is that it is possible to impose the condition that a 
is the only relevant Z 2 -odd operator, causing the allowed region in (A^., AJ-space to become 
a closed island, independent of auxiliary assumptions [9]. Because our A = 43 island lies 
within the error bars of [5], our results verify c-minimization to the precision achieved in 
[5].^^ Our results also give further evidence for the conjecture that the 3d Ising CFT is 
the unique Z 2 -symmetric 3d CFT with two relevant scalars and A^- < 0.6. (The precise 
condition on A^-, A^ depends on the shape of the allowed region further away from the 3d 
Ising point.) It would be very interesting to prove these conjectures analytically, perhaps 
by showing that the island in hgure 2 shrinks to a point as A —)■ oo. An alternative is that 
the conjectures are still true, but one needs information from other four-point functions to 
prove them. 

The analysis of [9] did not use permutation symmetry of the OPE coefficient = 
Ao-ecr-^^ Including this constraint leads to an additional modest reduction in the allowed 
region,which we plot in hgure 2 at A = 43. The resulting island gives a rigorous 
determination Ao- = 0.518151(6), A^ = 1.41264(6), which is 5-10 times more precise than 
the Monte Carlo results of [67]. We summarize the comparison to Monte Carlo in hgure 3. 


4 Discussion 


With SDPB, we have signihcantly improved the precision of (A^., A^) in the 3d Ising CFT. Our 
numerics indicate that the window of allowed dimensions may shrink to a point in the limit 
of inhnite computer time. In other words, they suggest that consistency of the correlators 

i2por readers interested in numerology, we recommend [68]. However, see [69]. 

^■^To see this more explicitly, it should be possible to place both upper and lower bounds on c using SDPB 
and see that it is constrained to lie close to the minimum computed in [5]. We leave this to future work. 

^^Note that permutation symmetry holds only when the conformal blocks are correctly normalized as 
functions of A. For scalars, the correct normalization is gA,o{u,v) = -I- ... to leading order in m, up to 

a A-independent constant. 

^®This fact was discovered during collaboration with Filip Kos, David Poland, and Alessandro Vichi [71]. 
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allowed region for A = 43, using three-point symmetry 


A, 

1.41285 : 
1.4128 : 
1.41275 : 
1.4127 : 
1.41265 : 
1.4126 : 
1.41255 : 
1.4125 L 



.-^^^^^^^^^^^^^^^^^ A,, 

0.51814 0.51815 0.51816 0.51817 


Figure 2: Allowed region for a Z 2 -symmetric 3d CFT with two relevant operators, computed 
with SDPB at A = 43. The light-blue region is a zoom of the smallest region in figure 1. The 
darker-blue region additionally uses symmetry of the OPE coefficients Xaae = ^aea- The black 
rectangle shows the estimate for (Aq-, A^) using c-minimization at A = 41 [5]. 


{{(TCTcrcr), (crcree), (eeee)}, together with the assumption that a and e are the only relevant 
scalars in the theory, may uniquely £x the dimensions (Ao-, A^). This conjecture could be 
more tractable analytically than trying to solve the full CFT consistency conditions. 

There are many more 3d Ising observables to explore. For example, the coefficient 
should be computable, and it would be interesting to compare with the recent Monte Carlo 
prediction [72]. It will also be important to consider larger systems of 3d Ising correlators. 

However, SDPB should also enable wider exploration of new correlators and diverse 
theories. SDPB is already being used in several bootstrap studies that would have previously 
been difficult or impossible [71, 73]. An exciting direction that may now be accessible is 
studying a four-point function of stress-tensors in 3d CFTs. 

In addition to the four-point function bootstrap, semidehnite programming has also 
recently been applied to the “modular bootstrap” in 2d CFTs [74, 75]. SDPB is equally 
applicable to modular bootstrap computations, since they too can be phrased in terms of 
polynomial matrix programs. 

From the computing point of view, there are many opportunities for improvement. For 
example, it should be possible to parallelize SDPB up to hundreds of cores, which could 
lead to even more precise calculations, and (just as importantly) easier exploration. Very 
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comparison to Monte Carlo 


A, 

1.4131 r 
1.413 : 
1.4129 : 
1.4128 : 
1.4127 : 
1.4126 : 
1.4125 : 
1.4124 ^ 


Monte Carlo 



L 


j_._I_._._._._^_._._._._^_._._._._^_._._._._I_._._._._^^ ^ 

0.51808 0.5181 0.51812 0.51814 0.51816 0.51818 


Figure 3: Comparison between the allowed region for the 3d Ising CFT using SDPB with 
A = 43 (blue) and Monte Carlo determinations of critical exponents (dashed rectangle) [67]. 
The size of the Monte Carlo rectangle is set by statistical and systematic errors associated 
with the simulation. By contrast, the blue region is a rigorous bound with sharp edges. 


different algorithms, like Second Order Conic Programming (SOCP), cutting plane methods, 
or constrained nonlinear optimization may also be applicable. 

The revival initiated in [1] is still young, and the technology (both analytical and 
numerical) is evolving rapidly. Current techniques are likely not maximally efficient, and it 
will be important to consider other methods, from new algorithms and optimization tools 
to conceptually different approaches. We are optimistic that much more will be possible. 
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A Choices and Parameters 


The PMP for {{aaaa), (aaee), (eeee)} in the 3d Ising CFT depends on the following choices: 


• An integer A specifying which derivatives to include in the functional a. These are 
given by with m> n and m + n < Ad® Because of the symmetry/antisymmetry 

of v) under u -H- n, some of these derivatives vanish identically. Including 

only non-vanishing derivatives, the dimension of a is 


dim a 




(A.l) 


• An integer k, controlling the accuracy of the approximation for conformal blocks (3.7). 
The positive prefactor is 


Xii^) 


-i- > 0 

n.(A-A„) - 


(A.2) 


where r* = 3 — 2^/2 is the radius of the point z = z = ^ in the radial coordinates of 
[76, 77]. The dimensions A** are special values below the unitarity bound, so that the 
product ni(A — A*j) is positive for unitary theories. The approximation (3.7) can 
be systematically improved by including more poles (A — A**)”^ and increasing the 
degree of pf‘^’^''*’™'"'(A). Our choice of poles is 


I 1 — i — k I k = 1,..., K, 

A*i e < d/2-k, I k = 1,... ,[k/2\, 

\^i + d —1 —k I /c = 1 ,..., min(K, [£/2j) 


(A.3) 


Smaller k means smaller-degree polynomials and shorter runtimes. Larger k is needed 
to get an accurate approximation for conformal blocks. We choose k by computing 
bounds with successively larger values of k until the results stabilize. Our hnal values 
are conservative: smaller k may still give sufficient accuracy. Derivatives of conformal 
blocks were computed using the recursion relation in [9] to order r®® (far greater 
accuracy than needed). 

• A set of spins S = {£i, ...,£/,} to include. If not enough spins are included, the solver 
may hnd a functional a that violates a positive semidehniteness constraint for some 
spin. Because derivatives of conformal blocks converge as a function of spin, in practice 
it is sufficient to include a hnite number of spins to ensure a satishes the constraints 
for all spins (as can be verihed post-hoc by testing a on constraints that were not 
explicitly included). This sufficient number of spins grows with A. Our choices are 
given in (A.4). 

^A=i9 = {0,..., 26} U {49, 50} 

Sa=27 = {0,..., 26} U {29, 30, 33, 34, 37, 38,41,42,45,46,49, 50} 

,Sa=35 = {0,..., 44} U {47,48, 51, 52, 55, 56, 59, 60, 63, 64, 67, 68} 

^^^43 = {0,..., 64} U {67, 68, 71, 72, 75, 76, 79, 80, 83, 84, 87, 88}. (A.4) 

= 2nmax — 1, where nmax is the parameter defined in [9]. 
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Once these quantities are fixed, the parameters to SDPB must be chosen to ensure numerical 
stability, precision, and correctness. Our choices for the computations in this work are 
summarized in table 1. 


A 

19 

27 

35 

43 

K 

14 

20 

30 

40 

spins 

*S'a=19 

Sa=27 

*S'a=35 

*S'a=43 

precision 

448 

576 

768 

960 

findPrimalFeasible 

True 

True 

True 

True 

findDualFeasible 

True 

True 

True 

True 

detectPrimalFeasibleJump 

True 

True 

True 

True 

detectDualFeasibleJump 

True 

True 

True 

True 

dualityGapThreshold 

0 

1 

CO 

0 

0 

1 

CO 

0 

0 

1 

CO 

0 

10-75 

primalErrorThreshold 

0 

1 

CO 

0 

0 

1 

CO 

0 

0 

1 

0 

10-75 

dualErrorThreshold 

0 

1 

CO 

0 

0 

1 

CO 

0 

0 

1 

0 

10-75 

initialMatrixScalePrimal {fl-p) 

1040 

1050 

1050 

IO6O 

initialMatrixScaleDual (flu) 

1040 

1050 

1050 

IO6O 

feasibleCenteringParameter (/^feasible) 

0.1 

0.1 

0.1 

0.1 

infeasibleCenteringParameter (Unfeasible) 

0.3 

0.3 

0.3 

0.3 

stepLengthReduction ( 7 ) 

0.7 

0.7 

0.7 

0.7 

choleskyStabilizeThreshold (6) 

10-40 

0 

1 

0 

lO-ioo 

10-140 

maxC omp1ement arity 

lOioo 

10130 

IOI60 

10200 


Table 1: Parameters for the computations in this work. Only SDPB parameters that affect the 
numerics (as opposed to parameters like maxThreads and maxRuntime) are included. The sets 
of spins 5 a are given in (A.4). Variables in the interior point algorithm of section 2.4.4 that 
correspond to SDPB parameters are indicated in parentheses, precision is in binary digits. 
The spin sets S\ refer to (A.4). 


B Performance 

B.l Complexity Comparison 

Let us compare the complexity of SDPB’s algorithm to that of SDPA-GMP for solving PMPs. 
For simplicity, suppose that each polynomial matrix M^{x) has size m x m and degree d. 
The matrices X, Y then have 2 J blocks, each of dimension m{d + 1). We focus on the most 
expensive parts of each iteration of the interior point algorithm and count the number of 
multiplications to leading order in J, m, and d. 

For SDPA-GMP, we assume the setup described in [17, 9]. There, the free variables y 
were embedded as 1 x 1 diagonal blocks in the matrix Y. (By contrast, for SDPB they are 
treated separately.) The most important contributions to the running time of SDPA-GMP are 
as follows. 
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• The Schur complement matrix is dense, so each of its elements must be computed indi¬ 
vidually. Computing X~^AqY requires N + 2J{mdY multiplications, since it involves 
a block-diagonal dense matrix multiplication X~^ x {AqY) (the linear term in N comes 
from multiplying 1x1 blocks). This must be repeated times: once for each q. 
Now, the matrices Ap are typically sparse, with 0{md) entries. Thus, evaluating all 
the traces TT{ApX~^AqY) requires approximately {N + 2J{md)^)^^^^+ 
multiplications. These steps dominate the running time for the computations in [9]. 

• Because the Schur complement matrix is dense, it must be inverted using a full 
Cholesky decomposition, which requires 1( -^^ multiplications. 

For SDPB, computing S takes negligible time. The most important steps are in solving 
the Schur complement equation: 

• Computing the Cholesky decomposition S = LL^ takes multiplications, 

since it can be done block-wise. 

• Forming L~^B requires N multiplications. 

• Forming Q = B) requires multiplications. This step domi¬ 

nates the running time for the computations in this work. 

• Computing the LU decomposition of Q requires multiplications. 

B.2 Running Time for 3d Ising Computations 

For the 3d Ising computations in this work, we have 

5 

d ~ A + -K (for large enough £), 
m = 1 or 2, 

J K, number of included spins. (B.l) 

A full analysis of the complexity in A would require determining the correct asymptotics 
of each quantity (including SDPB parameters like precision), which may depend on the 
computation at hand. In table 2, we simply report average runtimes and approximate 
multiplications/iteration for the choices given in appendix A. 
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solver 

A 

runtime (dual feasible) 

runtime (dual infeasible) 

mul./iter. 

SDPB 

19 

3.5 

0.9 

2 X 10« 

SDPB 

27 

32 

7.6 

1 X 10^ 

SDPB 

35 

190 

40 

5 X 10^ 

SDPB 

43 

810 

260 

2 X 10^° 

SDPA-GMP 

19 

~ 300 

~ 300 

1 X 10^^ 

SDPA-GMP 

27 

- 

■- 

1 X 10^2 

SDPA-GMP 

35 

~ 

- 

8 X 10^2 

SDPA-GMP 

43 

- 

- 

5 X 10^3 


Table 2: Runtimes for a single feasibility computation in the 3d Ising CFT using the 
correlators {(crcrcrcr), (crcree), (eeee)}, as described in [9] and appendix A. Average runtimes 
are different depending on whether a spectrum is disallowed (dual feasible) or allowed (dual 
infeasible). All times are in CPU-hours (for SDPB, this means the actual runtime is multiplied 
by maxThreads, which was 16 for most of the computations in this work). Approximate 
SDPA-GMP times are from [9]. All computations were performed on 3.3GHz 64-bit Intel 
Xeon Sandy Bridge processors. The column “mul./iter.” gives the approximate number of 
multiplications per iteration, calculated according to the discussion in subsection B.l. (To 
estimate running time from the number of multiplications per iteration, one needs to take into 
account precision, which also increases with A.) The SDPA-GMP computations with A > 19 
have not been attempted. 
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